Dynamical mean field solution of the Bose-Hubbard model 
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We present the effective action and self-consistency equations for the bosonic dynamical mean field 
(B-DMFT) approximation to the bosonic Hubbard model and show that it provides remarkably 
accurate phase diagrams and correlation functions. To solve the bosonic dynamical mean field 
equations we use a continuous-time Monte Carlo method for bosonic impurity models based on 
a diagrammatic expansion in the hybridization and condensate coupling. This method is readily 
f^*) , generalized to bosonic mixtures, spinful bosons, and Bose-Fermi mixtures. 
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Dynamical mean field theory (DMFT) is a computationally tractable framework for the study of fermionic lattice 
>> \ | models, which becomes exact in the limit of infinite dimensions or infinite coordination number [K-Q. In this limit 
the self-energy is momentum independent and can be obtained from the solution of an appropriately defined impu- 
rity problem. In finite dimensions, the approximation of a momentum-independent self-energy means that spatial 
correlations are neglected, but the local dynamics can be fully taken into account. DMFT provides a powerful, non- 
pcrturbative tool to investigate correlation effects and has been used extensively to study the properties of strongly 
correlated electron systems 0, 0] ■ The DMFT formalism is particularly well-suited to study the Mott transition, 
which is driven by local physics. 

The formulation of a dynamical mean field theory for bosonic lattice models is related to the formulation of 
an extended-DMFT for the ordered phase Attempts to derive bosonic DMFT equations based on the infinite 
'"^J . coordination limit face the problem that normal and condensed bosons must be scaled differently, rendering the 
£h ' scaling ambiguous (in fact the pathological infinite coordination limit can only be defined for a classical field Q). In a 
recent paper Byczuk and Vollhardt suggested to perform the scaling in the action rather than in the Hamiltonian. 
They obtained a conventional DMFT description with a momentum independent self-energy which requires the self- 
consistent solution of an impurity problem which couples to two baths, a bath of normal bosons and a bath of condensed 
bosons. While Ref. @ presented results for the bosonic Falicov-Kimball model, other groups applied (variants of) the 
B-DMFT formalism to the Bose-Hubbard model @ and multi-component Bose gases However, these previous 
formulations of B-DMFT were either incorrect or incomplete. 
. Here we derive internally consistent B-DMFT equations for the Bose-Hubbard model, which produce physically 
meaningful solutions over the whole parameter range, and recover the limits of the dilute Bose gas (loj and non- 
t^J- • interacting bosons. By comparing to the numerically exact lattice QMC simulation [ll[ of the full lattice model we 
show that our formalism yields remarkably accurate results [l2| both for phase diagrams and correlation functions. 
The DMFT equations for fermionic models can be derived (without any reference to an infinite coordination limit) 



using functionals of the local Green's function j 1 31 ] - with DMFT corresponding to a certain approximation of the 
kinetic energy functional. Our derivation of the B-DMFT equations is completely analogous [13]. 

In both the fermionic and bosonic versions of DMFT, the computationally challenging part is the solution of a 
quantum impurity problem. For fermionic impurity models, significant progress has been made with the development 
of continuous-time Monte Carlo techniques, based on an expansion of the partition function in powers of the interaction 
15 , 1(| or the impurity-bath hybridization [llj ■ In this Letter we show that a similar hybridization-expansion approach 



yields an efficient quantum Monte Carlo algorithm for bosonic impurity models. 

Our specific model is that of spinless bosons on a three-dimensional (3D) simple cubic lattice with Hamiltonian 

H = -t J2 bib, + |l>(n i -l)- j uX>, (1) 

where t denotes the hopping amplitude, U the on-site interaction and fj, the chemical potential. 

The effective impurity model of B-DMFT contains the microscopic local terms of the Hamiltonian (U and fx) 
to which two additional source fields are added. The first one can, like in static mean-field theory, be written as 
~zt(j> (with cj) the constant condensate) and is conjugate to the b and operators, such that (b) can become a non- 
zero complex number. Fluctuations around the condensate at the one-loop level can be added after decomposing 
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6(r) = (b) + Sb(r). The source field for the two-particle channel couples to non-condensed operators [19( of the form 
5b^ (T)Sb(T 1 ) (and also 5b(T)Sb(r') because of the symmetry breaking). This follows from a functional derivation [l4[ 
in which the connected local Green's function is constrained [l8| . If the full local Green's function were constrained, 
one would not arrive at the correct B-DMFT equations, because the approximation of the kinetic energy functional 
would mix up small and large contributions. 

Writing the partition function as Z = Tr;,[Te~' Simp ] and shifting contributions between the hybridization and 
condensate terms, our final action can be expressed using the full operators b in the (Nambu) form: 
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1 fP fP 

- dTdr'h^ (t)A(t — r')b(r') — M / d,Tn{r) 

2 Jo Jo 

/ dm(r)[n(r) - 1] - k& / drb(r). (2) 
Jo Jo 



The boson creation/annihilation operators are = (b' , b), the time- independent condensate is = (</>*, </>), and the 
hybridization function A is related to the mean-field propagator Go through 

A(iu n ) = -«w„tr 3 -p,l + Gq (3) 

The parameter /Et = ji — (e) is chosen such that A(zu;„) — > in the limit cj n — > oo. We will consider here a symmetric 
density of states, where (e) = 0, and just write fj, from now on. We furthermore define the elements of the hybridization 
matrix as 

*r-0-(,S?- T i^r-J?V (4) 



2K*(t-t') F(t-t') 

The condensate $ is constant in time and determined by the simple self-consistency condition 

* = (b(r)) Simp . (5) 

To determine the hybridization function A we calculate the matrix self energy using the Dyson equation 

S( iW „) = G- 1 (^„)-G- 1 ( iW „), (6) 

where the connected part of the Green's function for the normal bosons is given by G c (r) = — (Tb(r)b^(0)),s im +<!><I?t. 
From X and the dispersion of the lattice we obtain the local lattice Green's function 



Gl a tt(«W„) = ^ [ iuJ nO~3 + - - £(iw n ) , (7) 



and the self-consistency condition requires that the impurity Green's function coincides with the local lattice Green's 
function: 

Gq 1 ^) = -E(ioj n ) + G^ tt {iu n ). (8) 



Assuming K = K* and <f> = <j>* , the coupling k between the impurity and the condensate is given by 14 1 

K = zt - F(iu n = 0) - 2K{iuj n = 0). (9) 

If one incorrectly uses k — zt the trivial limit of the dilute Bose gas [l(| and the non-interacting Bose gas is 
incorrectly reproduced in finite-dimensional systems, while the condensed phase becomes unstable in a wide region of 
parameter space. Hu and Tong [f| set the off-diagonal hybridization function K to zero after each iteration to avoid 
these instabilities, and thus solved incorrect B-DMFT equations, which equally fail to reproduce the trivial limits. In 
Ref. B-DMFT was considered as an expansion in l/z on a tree-like structure and a perturbative prescription was 
used to rescale the condensate (or k) after each step. Complementing their l/z expansion approach with the correct 
non-perturbative prescription would lead to Eqs. and ^ f20j ] . 

Our action reproduces the correct results in all limits on any lattice, and also produces stable B-DMFT solutions 
including in phases with a non-zero condensate. In the non-interacting model, the chemical potential is pinned at the 
lower band edge in the presence of a finite condensate and Eq. © reduces to k = — Gq 1 (iui n = 0). Note that there 
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FIG. 1: Diagram corresponding to perturbation orders mr = 1, wik = 1, mic* = 1, "i<p = 2, m^,* = 2 and n(r = 0) = 2. 



is an isolated state at energy — zt for the Bethe lattice DOS [21[. In the static case without hybridization Eq. © 
gives k = zt consistent with static mean field theory. In the limit of infinite dimensions, only static mean-field theory 
is physical, in which case the addition of the source field A is impossible because 8b = 0. In fact, the entire phase 
diagram of the 3D Bosc-Hubbard model, as well as dynamical quantities such as correlation functions, are reproduced 
by the B-DMFT action with remarkable accuracy. 

The self-consistency equations are solved by starting from an initial guess, solving the quantum impurity problem 
and then calculating new values for <fr through Eq. ([5]), A through Eq. ([H]) and k through Eq. ©. This procedure is 
repeated until convergence is reached. 

The computationally demanding step is the solution of the bosonic quantum impurity problem and we will now 
present a quantum Monte Carlo (QMC) algorithm for its solution, which is similar in spirit to the fermionic hy- 
bridization expansion algorithm of Ref. (l7j . We expand the partition function Z = Tr/,[Te _ ' Simp ] in powers of the 
hybridization functions F, K , K* and the source fields <fi and 0*. This leads to an expression for the partition function 
as a sum of diagrams of the type illustrated in Fig. [TJ which can be represented by a collection of mp + Iran- + m^ 
creation operators and the same number mp + 2m k + of annihilation operators on the imaginary time interval 
[0,/3). Hybridization functions F connect mp pairs of creation and annihilation operators, off-diagonal hybridization 
function K (K*) connect mjf (rriK*) pairs of creation (annihilation) operators, while m<f, (m^*) creation (annihilation) 
operators are linked to source fields </> (0* ) . The integer n > corresponds to the occupation of the impurity at times 
t = and ft, and thus fixes n(r). 

An ergodic sampling of all possible diagrams requires the following updates: (i) insertion/removal of a pair b(r)F{T — 
T')b^(r'), (ii) increase/decrease of n by one, (iii) change of the bath type: 



b{r)F{T - T')b\r') O K^b^)^ (r'), 
b{r)K*{T - r')6(r') O k0*6(t)k0*6(t'), 
b\r)K{T - t'P{t') o K<t>b\T)K<t>tf{T'). 



(10) 

(11) 
(12) 



Additional updates such as shifts of operator times and reconnections of hybridization lines can be used to improve 



the efficiency. Denoting the trace contribution of a diagram (n\ . 
detailed balance condition for inserting/removing a pair b(r)F(r 



. \n) by WTr(n;T{ , 
T')tf(T') becomes 
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if we choose the times r and r' randomly in [0,/3) and propose to remove this pair with probability l/(mp + 1). 
Similar formulas can be derived for the other updates. 

By taking functional derivatives of the partition function with respect to cither F(t — r'), K(t — r'), K*(t — t'), 
or 0* and 4> one obtains measurement formulas for the diagonal and off-diagonal parts of the Green function matrix 
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FIG. 2: Top panel: phase diagram (superfluid to normal liquid transition) in the space of interaction and temperature for 
n = 1. The dashed line shows the static mean field result, the red curve the exact solution for a Bethe lattice with coordination 
number z = 6 (Ref. l2(il) and the blue curve with open diamonds the QMC result from lattice simulations (Ref. fl2l) . The black 
line with open circles corresponds to the B-DMFT solution, which yields a second order transition. Bottom panel: ground-state 
phase diagram in the space of t/U and fi/U, showing the first two Mott lobes surrounded by superfluid. The B-DMFT phase 
boundary was computed at /3t = 2. Error bars are much smaller than the symbol size. 



and the condensate order parameter: 



- (gM^rL (14) 

Wr)Kff)>«. - If. ■ («) 



MC 




(b(r))s itap = V . (16) 



and similarly for the adjoint with A(r, f ) = 5(t — f) for f > and 5(t — f — /3) for f < 0. 

The end point of G(r) = — (Tb^r)^ (0)) s imp is given by the density G(/3_) = — (n) MC = —(4 Jq drn(r))MC an d 
G(0 + ) = G(/3_) — 1. Here (A)mc means that the quantity A should be averaged over all configurations obtained 
in the Monte Carlo sampling. Although the functions K{t) are negative (well known for the weakly interacting 
Bose gas [13]) and lead to a sign problem, this only becomes an issue at very low temperatures in the presence of a 
condensate (see Tab. HJ but does not prevent an accurate computation of phase diagrams and dynamical quantities. 
For the temperatures considered here one iteration takes just a few minutes on a single CPU core. 

As main results we show the finite temperature phase diagram (top panel of Fig. [5]) and the ground state phase 
diagram (bottom panel of Fig. [2]) for the first and second lobe of the Bose Hubbard model on a 3D simple cubic lattice 
and compare results obtained with B-DMFT to exact results from lattice QMC simulations, the exact solution for the 
Bethe lattice with coordination number z = 6 [2fj| . and to static mean field results. For the calculation of the ground 
state phase diagram we used fit = 2, which is shown in Fig. [2] to be a sufficiently low temperature [2Sj] . 

The excellent agreement between our B-DMFT results and the full solution of the Bose-Hubbard model shows 
that the Mott-transition is a local phenomenon, well described by a momentum-independent self-energy and that 
the condensed bosons are accurately described by a uniform condensate. B-DMFT also yields remarkably accurate 
data in the condensed phase, as illustrated in Tab. HI Even dynamical quantities such as density-density correlation 
functions are correctly reproduced in all three phases (Fig. [3]). 

The generalization of our QMC algorithm to bosonic mixtures is straightforward. It can also be extended to spinful 
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FIG. 3: Comparison of the connected density-density correlation functions (n(r)n(O)} — (n) 2 in the normal ([7=20, /x=8.72, 
,3=0.2), superfluid ((7=10, /i=0.7, /3=2) and Mott insulating (17=40, /x=ll, /3=2) phases. The left panel shows the B-DMFT 
result and the right panel the exact data from lattice QMC. The data for the Mott phase are scaled by a factor 10. 



bosons with spin-dependent or more complicated interactions using a matrix-formulation analogous to the fermionic 
algorithm of Ref. [22|. Similarly, the method can be extended to impurity clusters to solve cluster-generalizations 23 1 of 
B-DMFT. These extensions will enable the simulation of superfluids, supersolids, and super-counter-fluids in bosonic 
mixtures and spinor condensates and insulators in spinful bosonic systems. 

Our bosonic impurity solver, based on a hybridization expansion, can easily be combined with its fermionic coun- 
terpart [it]], enabling an efficient DMFT simulation of Bose- Fermi mixtures. While powerful numerical methods exist 



for the simulation of bosonic lattice models [III ] , these simulations become prohibitively difficult as soon as fcrmions 
are involved. We thus view the solver and formalism described here not only as an efficient tool for the solution of 
bosonic problems, but also as an important step which opens the door to a systematic investigation of Bose-Fermi 
mixtures [24j |. 
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sign 


(20, 6.6) 


0.99441(4) 


0.99456(1) 


0.5042(3) 


0.486(2) 


0.6373(1) 


(24, 8.6) 


0.99494(5) 


0.995120(1) 


0.3383(4) 


0.316(1) 


0.7836(1) 


(26, 10) 


1.00194(3) 


1.001936(1) 


0.2389(4) 


0.2227(9) 


0.8674(1) 


(28, 11.3) 


1.00252(3) 


1.002598(4) 


0.1087(5) 


0.104(1) 


0.9585(1) 


(30, 13) 


1.000403(5) 


1.00041(4) 








1 


(32, 15) 


1.000333(5) 


1.000370(9) 
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TABLE I: Comparison of n and 4> 2 between lattice QMC and B-DMFT for different values of /i and U at fit = 1. 
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Auxiliary material for Dynamical mean field solution of the Bose-Hubbard model 



Scope 



In this supplementary material we provide a formal basis for the action and the self-consistency relations that were 
presented in the manuscript. We implement an expansion around the atomic limit, following almost literally the 
lecture notes by A. Georges [l3j . and consider B-DMFT as an approximation to the kinetic energy functional. The 
atomic reference system is interpreted as our impurity problem. Use is made of the coupling constant integration 
method and source fields (Lagrange multipliers) are introduced to constrain the condensate field and the connected 
Green's function for the normal bosons to their physical values. 

Expansion parameter 

We introduce a parameter a €E [0, 1] such that denotes nearest neighbor sites) 

When a = 0, the atomic limit is recovered and the partition function factorizes over all sites. When a = 1 the full 
hopping is recovered, and this is the model we are ultimately interested in. 

Source fields and constraining fields 

Constraining the normal/anomalous Green's functions and the condensate to specified values can be done by 
introducing conjugate source fields (Lagrange multipliers) in the action. In order to constrain the condensate to <& we 
introduce the source field J, and analogous for the connected Green's function G c with source field A. Throughout 
this document we use the Nambu notation in which = (</>*, 0), = (J*, J) and the individual components of G c 
and A are given by 

and 

*m - ( F i~:\ 2K{ 1 ) w 

V ' y 2K*(t) F(t) ) v 1 

We can then explicitly write down the grand potential per site (there are N s sites) which is a functional of the source 
fields and also depends on the constraining fields, 



r! Q [J,#,A,G c ] = 



l — In J V[b*,b] exp | J* dr (E K(-dr + n)h - H a [b*,b]j 







+ / drE(^W^(-r)-^(r)] + J(r)[6|(r)-^(r)]) 
Jo 



l dT f dT'y2F(T~T')[5h(T)5b*(T')+G c (T-T')} 

Jo Jo i 

f dT f dT' E K(t - T ')[6b*(T)5b*(T') +G* c (t- r')] 
Jo Jo l 

j* dT J* dT 1 E K*(t - r')[^(r)^(r') + G c (t - t')]\. 



(20) 
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Atomic limit : impurity model 

Let us consider the case a = 0, which is the atomic limit. The problem becomes local on every site with grand 
potential 



r! [J,*,A,G c ] = - 



l — In J V[b*,b] cxp | J* dr Ki-dr + n)h - ^(n, - 1)J 



drJ2 (Jo (r)[h(r) Mr)} + Jb(r)[6J(r) - (t)]) 

i 

f dr f dr' VF„(r- r'^Sb^SbUr') + G c (r- r')\ 
Jo Jo i 

f dr f dr' J2 Ko(r r')[5b*{r)8b*{r') + G* c (r - r')] 
Jo Jo i 

J* dr J* dr 1 J2 K* (r - T')[Sb t (r)Sb t (r') + G c {r - r')]|. 



(21) 



From SQq/SJq = and 8Qo/5Jq = we obtain 

* = (b)s lmp , (22) 

and from 5£l /dF = 0, 5£l /5K = 0, SCl /SK* = the relation 

G c (r) = G(r) + **t : (23) 

with G(r) = -(Tb(r)bt(0)) Simp . The expectation values (. . .) Simp = Tr[Te~ s ^ . . .}/Tr[T e - s ^} are defined with 
respect to an impurity action 



\ f f drdr'5b*(r)A (r-r')6b{r')- f drjj (r)b(r) 
2 Jo Jo Jo 



"Jimp — 2 



Jo Jo 
rP u f P 

-fx drn(r) + — drn(r)[n(r) - 1]. (24) 
Jo 2 7 

Inverting expressions (|2"2"|) and ((23)) yields J[$,G C ] and A [3?,G C ] and thus a functional T of the condensate and 
connected impurity Green's function: 

ro[*,G c ]=F imp [*,G c ]-^dT[Fo(r)G c (r) + Xo(r)G:(r) + A *(T)G c (r)] + ^ / V^[J *(t)<Mt) + Jo(r)^)]. 

(25) 

Full model 

The exact functional of the (local) Green's function and condensate are constructed using the coupling constant 
integration method, starting from the atomic limit, T = T a —i — To + j da^-^-. Using the stationarity of Q, (a- 
derivatives of the Lagrange multipliers do not contribute), 

rfp i rP i rP 

~d = ~JTB / dTt T,^ T ) h ^ T )) = -JTs / dT *S^*( r )^( T ) + ( 6b *(r)Sbj(r))} (26) 

B r J (J i ■ - \ S J (J / - - \ 

= Fs^E^cCk,^)!^--^ f drtJ2^t(r)Mr)}- (27) 

We arrive at the formal expression for the exact functional T = r Q= i, 

r[*,G c ] =r [*,G c ]+/C[*,G c ], (28) 
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with the kinetic energy functional /C[3>, G c ] = J da^f [3?, G c ]. Requiring stationarity (5T/84>*(t) = 0, <5r/<5 </>.,■ (t) = 
0) determines the value of the source field conjugate to the condensate (assume a homogeneous condensate over the 
lattice), 

J = zt<f>. (29) 

Since the condensate is time- independent (and taken real), we drop the r dependence of Jo as well. The other 
stationarity requirement (ST/6G C = 0, 5T/SG C = 0) determines the hybridization function: 

Fo(r) = J^- y K Q {r) = -^- (30) 
5G c (t) 5G c (t) 

Note that for the case z = oo, we have identically 8b = and only static mean- field theory exists (fjj. 

Approximation to the kinetic energy functional 

B-DMFT can now be considered as an approximation to the kinetic energy functional. With the single-particle 
Green's function of the Bose-Hubbard model in the presence of source fields and for arbitrary coupling constants, we 
can define a sclfencrgy 

G^(k,iu> n ) = [iu n a 3 + (/i - aek)I + A Q [iw n ] - S Q [k, iw n ]] _1 . (31) 

The DMFT approximation consists in replacing the self-energy S a for arbitrary a by the impurity model self-energy 
S - Hence, 

G"(k, iw„) | b-dmft = [iw„cr 3 + (// - ae k )I + A Q [iw„; *, G c ] - S Q= oN«; *,G C ]] _1 (32) 
= [A a - Ao + G-i-aetl]- 1 , (33) 

where we have used that the impurity self-energy satisfies the Dyson equation S Q= o[iw„; G c ] = iLo n <r^ + /xl + 
Ao[iw„; G c ] — G~ . Summing over k, and using the constraint on the local lattice Green's function, we obtain the 
following relation between G c and the hybridization function: 

G c (*<) = [ deD{e)(( - ael)- 1 = -d(-), (34) 
J a \a/ 

with ( = A a — A + Gj 1 . We used the non-interacting density of states is D(e) — J2 k 8(e — £k) and its Hilbcrt 

transform D(z) = J deD(e)(z — el) . By introducing its inverse, D(R(g)) = g, the relation above can be inverted 
(aR(aG c ) = ( = A Q — Ao + G" 1 ) and yields the hybridization function as a functional of the local Green's function, 

A a [icj n ; G c ] = -G" 1 + A [*, G c ] + aR(aG c ). (35) 

The lattice Green's function expressed as a functional of G c becomes 

G"(k, i^-Vi) | b-dmft = (aR(aG c ) - aekl)^ 1 - (36) 

Equation (|27|) can now be evaluated with G" (k) | b-dmft : 

rP 

-Tr f,,G a Ck iw_ W-o t~\~\ h?t — I 



L T r^e k G«(k ) io; n )|B-DMFT-^-g / drt $>?(r)&(r)] 
n,k s 

= -^E Tr / deeD(e)(R(aG c )-eI)-'--±- f dr^[^(r)^(r)]. 
= \\ E ^ [ - 1 + «G c i?(aG c )] " ^ drt E ( T )^ ( T )l ' ( 37 ) 



An explicit expression for the B-DMFT approximation to /C[<&, G c ] therefore reads 



/C B -dmft[*,G c ] = / dcJ-^Tr [G c {iu n )R{a.G c {iUn)) ~ a" 1 !] - J- ( drt ^[^(r)^(r)], (38) 



where the last term reduces to —zt<p*<f) for a constant, homogeneous condensate. 
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Stationarity conditions 

It immediately follows from Eq. (f38|) that the stationarity condition for the condensate is unaltered in the B- 
DMFT approximation (Jo = zt&), while the stationarity condition for the connected Green's function (ST/SG C = 0, 
ST/5G C = 0) reads in the B-DMFT approximation (use R(aG) + aGR'(aG) = d a [aR(aG)] and the cyclical properties 
of the trace), 

A [iw„;$,G c ]|b-dmft = -R[G c (ioj n )] +G c (iu n y 1 = -icj n a 3 - /il + S. imp + G,? 1 . (39) 
Applying £)(.) to both sides of Eq. ([5^)1 gives 

G c {iuj n ) = J deD(e)(iuj n a 3 + (fj, - e)I - E irap (za; n )) _:L . (40) 

This equation defines the B-DMFT self-consistency condition. 



Summary 

B-DMFT maps the bosonic lattice problem to a self-consistent solution of an impurity model, whose action (ex- 
pressed in terms of the full operators b) reads 



&-dmft = -■= / / dTdT'b*(T)Ao(T-r')b(r')-M / drn(T) + - d™(r)[n(r) - 1] 
^ Jo Jo Jo * Jo 

& (zt - dr'Ao(r')) ^ drb(r). (41) 

Using the fact that K = K* can be chosen real, one recovers the action presented in the manuscript. The solution 
of the impurity problem yields the condensate €> (Eq. (|22)l). the connected Green's function G c (Eq. (|23j)) and the 
self-energy Si mp of the impurity model. The right hand side of Eq. (|4"0")) then defines the local lattice Green's function, 
which is identified with the impurity Green's function and thus allows to define the new hybridization function for 
the next iteration using Eq. (|39j). 
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